Suppose \(Y_t\) is AR(1) series observation
(unconditional) Mean \[ E( Y_t ) = 0 \]
(unconditional) Variance \[ V( Y_t ) = \gamma(0) = (1+\phi_1^2) \sigma^2 \]
Conditonal Mean:
\[
E\Big(Y_t \hspace2mm \Big| \hspace2mm \mbox{all variables realized by yesterday. }\Big)
\]
Conditonal mean of AR(1): $Y_t = Y_{t-1} + e_t $ \[ \begin{align} E\Big(Y_t \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big) &= E \Big(\phi Y_{t-1} + e_t \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big)\\ \\ &= \phi Y_{t-1} + E(e_t) \hspace2mm = \hspace2mm \ph Y_{t-1}. \end{align} \]
Conditonal variance \[ \begin{align} V\Big(Y_t \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big) &= V\Big(\phi Y_{t-1} + e_t \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big)\\ \\ &= V(e_t) \hspace2mm = \hspace2mm \sigma^2 \end{align} \]
n=200
Y <- arima.sim(n = n, list(ar = c(0.8) ))
v.AR1 = (1+.8^2)*1
plot(Y, type="o", main="AR(1)", xlim=c(0,n*1.1))
abline(h=0)
abline(h=c(1,-1)*1.96*v.AR1, col="red")
lines(n+1, .8*Y[n], type="p", col="blue")
lines(c(n+1,n+1), .8*Y[n]+c(1.96,-1.96), type="p", col="red")Suppose \(Y_t\) is MA(1) series observation
(unconditional) Mean \[ E( Y_t ) = 0 \]
(unconditional) Variance \[ V( Y_t ) = \gamma(0) = (1+\theta_1^2) \sigma^2 \]
Conditonal mean of MA(1): \(Y_t = e_t + \theta_1 e_{t-1}\) \[ \begin{align} E\Big(Y_t \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big) &= E \Big( e_t + \theta_1 e_{t-1} \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big) \\\\ &= E ( e_t )+ \theta_1 e_{t-1} \hspace2mm = \hspace2mm \theta_1 e_{t-1} \end{align} \] Note that \(e_{t-1}\) is not observable.
Conditonal variance of MA(1): \(Y_t = e_t + \theta_1 e_{t-1}\) \[ \begin{align} Var\Big(Y_t \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big) \ &= Var\Big( e_t + \theta_1 e_{t-1} \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big)\\ &= Var( e_t ) \hspace2mm = \hspace2mm \sigma^2. \end{align} \]
AR(1) \[ \mbox{ Uncond'l } \hspace20mm \mbox{ cond'l } \\ E(Y_t) = 0, \hspace28mm E(Y_t|\omega_{t-1}) = \phi_1 Y_{t-1}, \\ Var(Y_t) = (1+\phi_1^2)\sigma^2 \hs10mm Var(Y_t|\omega_{t-1}) = \sigma^2 \]
MA(1) \[ \mbox{ Uncond'l } \hspace20mm \mbox{ cond'l } \\ E(Y_t) = 0, \hspace28mm E(Y_t|\omega_{t-1}) = \theta_1 e_{t-1}, \\ Var(Y_t) = (1+\th_1^2)\sigma^2 \hs10mm Var(Y_t|\omega_{t-1}) = \sigma^2 \]
For ARMA(p,q) model, conditional mean changes, but conditional variance is constant.
\[ \begin{align} X_t &= \mbox{ Stock Price (observation) } \\ \\ Y_t &= \ln(X_t) - \ln(X_{t-1}) \hs5mm : \mbox{ log return } \end{align} \]
library(quantmod)
source("http://gozips.uakron.edu/~nmimoto/477/TS_R-90.txt")
getSymbols("AAPL") #- download from Yahoo!## [1] "AAPL"
## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series
## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0.07 0.006 0.008 0 0 0 0.028
## [1] "GSPC"
## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series
## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0.001 0.001 0.003 0 0 0 0.007
Uncorrelated
Squares are correlated
Clustering
Asymmetry
Heavy Tailed unconditional and conditional distribution
Don’t confuse the conditional heteroscedasticity with (unconditonal) heteroscedasticity:
D <- read.csv("http://gozips.uakron.edu/~nmimoto/pages/datasets/gas.csv")
D1 <- ts(D[,2], start=c(1956, 1), freq=12)
plot(D1, type='o')Engle (1985) AutoRegressive Conditionally Heteroscedastic Model \[ \begin{align} Y_t &= \sigma_t e_t \hs10mm e_t \sim_{iid} \cN(0,1) \\\\ \sigma_t^2 &= \omega+ \alpha Y_{t-1}^2 \end{align} \]
Mean \[ E(Y_t) = \sigma_t E(e_t) = 0 \]
Variance \[ V(Y_t) = V(\sigma_t)V(e_t) = \]
Conditional mean \[ E[Y_t \Big| Y_{t-1}, e_{t-1}, \ldots] \hspace3mm = \hspace3mm \sigma_t E[e_t] = 0 \]
Conditional variance \[ V[Y_t \Big| Y_{t-1}, e_{t-1}, \ldots] \hspace3mm = \hspace3mm \sigma_t^2 V[e_t] = \sigma_t^2 \]
Conditional Mean: \(0\)
Conditional Variance : \(\sigma_t^2\)
Daily Price of SP500 ETF (SPY) from Jan 02 2000 to Dec 31 2014
D <- read.csv("http://gozips.uakron.edu/~nmimoto/pages/datasets/SPY.csv", header=T)
source("http://gozips.uakron.edu/~nmimoto/477/TS_R-90.txt")
head(D)## ixDay Date Weekday X SPY.Open SPY.High SPY.Low SPY.Close
## 1 1 1/3/2000 Monday 1 148.25 148.25 143.88 145.44
## 2 2 1/4/2000 Tuesday 2 143.53 144.06 139.64 139.75
## 3 3 1/5/2000 Wednesday 3 139.94 141.53 137.25 140.00
## 4 4 1/6/2000 Thursday 4 139.62 141.50 137.75 137.75
## 5 5 1/7/2000 Friday 5 140.31 145.75 140.06 145.75
## 6 6 1/10/2000 Monday 6 146.25 146.91 145.03 146.25
## SPY.Volume SPY.Adjusted
## 1 8164300 110.33
## 2 8089800 106.01
## 3 12177900 106.20
## 4 6227200 104.50
## 5 8066500 110.57
## 6 5741700 110.94
layout(1)
# Fit GARCH model
library(fGarch)
Fit1 <- garchFit(~ garch(1,1), data=X1, cond.dist="norm", include.mean = FALSE, trace = FALSE)
## Fit1@fit$par #-- estimated parameters
## Fit1@residuals #-- this is not the garch residuals!!!! same as X1
## Fit1@sigma.t #-- estimated sig_t
## Fit1@residuals/Fit1@sigma.t #-- this is the (standardized) GARCH residuals
print(Fit1@fit$ics) #-- AIC and BIC are here## AIC BIC SIC HQIC
## -6.470319 -6.465360 -6.470321 -6.468556
plot(X1)
lines( 1.96*Fit1@sigma.t, type="l", col="red")
lines(-1.96*Fit1@sigma.t, type="l", col="red")## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series
## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0 0 0 0.719 0.855 0 1.001
\[ \begin{align} \sigma_t^2 &= \omega + \alpha Y_{t-1}^2 + \beta \sigma^2_{t-1} \\\\ &= \omega + \alpha Y_{t-1}^2 + \beta \Big(\omega + \alpha Y_{t-2}^2 + \beta \sigma^2_{t-2} \Big)\\ \\ &= \omega + \beta \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \sigma^2_{t-2} \\ \\ &= \omega + \beta \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \Big(\omega + \alpha Y_{t-3}^2 + \beta \sigma^2_{t-3} \Big) \\ \\ &= \omega + \beta \omega + \beta^2 \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \alpha Y_{t-3}^2 + \beta^3 \sigma^2_{t-3} \end{align} \]
Continuing, we get \[ \begin{align} &= \omega (1 + \beta + \beta^2 + \cdots ) + \alpha \sum_{i=0}^k \beta^i Y_{t-1-i}^2 + \beta^{k+1} \sigma^2_{t-1-k} \\\\ &= \frac{\omega}{1-\beta} + \alpha \sum_{i=0}^\infty \beta^i Y_{t-1-i}^2 \end{align} \]
\[ \sigma_t^2 \hspace3mm = \hspace3mm \frac{\omega}{1-\beta} + \alpha \sum_{i=0}^\infty \beta^i Y_{t-1-i}^2 \] We will use the truncated and estimated version of this, \[ \hat \sigma_t^2 \hspace3mm = \hspace3mm \frac{\hat \omega}{1- \hat \beta} + \hat \alpha \sum_{i=0}^{t-1} \hat \beta^i Y_{t-1-i}^2 \]
\[ Y_t = \sigma_t e_t \]
Using observation \(\{Y_1, \ldots, Y_n\}\), the residuals are \[ \hat e_t = Y_t/ \hat \sigma_t \]
\[ \begin{align} Y_t &= \sigma_t e_t \\\\ \sigma_t^2 &= \omega + \alpha Y_{t-1}^2 + \beta \sigma^2_{t-1} \end{align} \]
[Conditonal distribution] \(=\) [Distribution of \(e_t\)]
Guessing correct distribution for \(e_t\) is very important in GARCH parameter estimation.
x <- seq(-10,10, .1)
#- Standardized t-distribution
plot(x, dstd(x, mean = 0, sd = 1, nu = 5), type="l" )
lines(x, dnorm(x, mean = 0, sd = 1), col="red" ) #- Generalized Error Distribution
plot(x, dged(x, mean = 0, sd = 1, nu = 4), type="l" )
lines(x, dnorm(x, mean = 0, sd=1) , col="red" ) #- Skewed-Standardized t-distribution
plot( x, dsstd(x, mean = 0, sd = 1, nu = 5, xi = 1.5), type="l" )
lines(x, dnorm(x, mean = 0, sd = 1), col="red" ) #- Skewed-Generalized Error Distribution
plot(x, dsged(x, mean = 0, sd = 1, nu = 5, xi = 1.5) , type="l" )
lines(x, dnorm(x, mean = 0, sd = 1), col="red" )#--- 5. Simulation and estimation of GARCH -------------------------------------------------
n <- 1000
theta <- c(.024, .1, .8) #- Parametes of GARCH(1,1)
#- Specify parameters and conditional distribuiton (distribution of e_t)
spec <- garchSpec(model = list(omega=theta[1], alpha=theta[2], beta=theta[3]), cond.dist="norm")
#- Simulate GARCH using above spec
x <- garchSim(spec, n = n, extended=FALSE)
#- Estimate parameters of GARCH
Fit1 <- garchFit(~ garch(1,1), data=x, cond.dist="norm", include.mean = FALSE, trace = FALSE)
coef(Fit1) #- estimated parameters
#--- Monte Carlo Estimation to see parameter MSE --
n <- 1000
theta <- c(.1, .1, .8) #- Parametes of GARCH(1,1)
#- Specify parameters and conditional distribuiton (distribution of e_t)
spec <- garchSpec(model = list(omega=theta[1], alpha=theta[2], beta=theta[3]), cond.dist="norm")
Param <- matrix(0, 1000, 3)
for (i in 1:1000) {
print(i)
Y <- arima.sim(
x <- garchSim(spec, n = n, extended=FALSE)
Fit1 <- garchFit(~ garch(1,1), data=x, cond.dist="norm", include.mean = FALSE, trace = FALSE)
Param[i,] <- coef(Fit1) #- estimated parameters
}
SE <- (Param - t( matrix(theta, 3, 1000) ))^2
MSE <- apply(SE, 2, mean) #- MSE of each parameter
rMSE <- sqrt(MSE); rMSETrue parameter (.024, .1, .8)
True cond’l dist: Normal Estimated using: Normal
True cond’l dist: std(5) Estimated using: Normal
True cond’l dist: sged(skew=.7, shape=1.45) Estimated using: Normal